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Abstract 

We compare nonlinear stresses and temperatures for adiabatic shear flows, using up to 262,144 
particles, with those from corresponding homogeneous and inhomogeneous flows. Two varieties 
of kinetic temperature tensors are compared to the configurational temperatures. This compar- 
ison of temperatures led us to two new findings, beyond our original goal of analyzing shear 
algorithms. First, we found an improved form for local instantaneous velocity fluctuations, 
as calculated with smooth-particle weighting functions. Second, we came upon the previously 
unrecognized contribution of rotation to the configurational temperature. 
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I. INTRODUCTION 



Hoover, Hoover, and Petravioi studied the nonlinear stresses and temperature changes 
induced by shear in a variety of stationary flows. They used nonequilibrium molecular 
dynamics to compare several popular algorithms. The algorithms are described briefly in 
Sec. H. See Figs. 1 and 2 for the geometries used to induce the flows. This work used a 
smooth repulsive soft-sphere potential^ with a range of unity, 

0(r < 1) = 100(1 -r^)^ . 

The particle mass, energy per particle, and density were all chosen equal to unity. These 
conditions correspond to a dense fluid at about 2/3 the freezing pressure. 

m = 1 ; E/N = {K + ^)/N = 1 ] p = Nm/V = 1 . 

Thermostat or ergostat forces were used to generate stationary states. There, as well as in 
the present work, we choose x for the flow direction and e = 0.5 for the strainrate, where 
the time-averaged velocity component of the flow Vx increases linearly in the y direction: 

Three-dimensional homogeneous periodic simulations [Fig. 1 shows a two-dimensional 
version] gave 

Tyy > T^x > T^z ; Pyy > Pxx > Pzz [Doll's Algorithm] ; 

Txx > Tyy > ; Pxx > Pyy > Pzz [SUod Algorithm] . 

and differed qualitatively from corresponding three-dimensional boundary-driven results 
[Fig. 2 shows a two-dimensional version]: 

Txx > T^z > Tyy ; Pxx > Pzz > Pyy [Bouudary - Driven] . 

The main conclusion drawn from that work was that neither homogeneous method, 
Doll's- nor Sllod4, successfully reproduces the more-physical boundary-driven results^. To 
quote Ref. 1: "The Doll's and Shod algorithms predict opposite signs for this normal- 
stress difference [Pxx — Pyy], with the Shod approach deflnitely wrong, but somewhat 
closer to the (boundary-driven) truth." 
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Figure 1: Two-dimensional version of periodic homogeneous isoenergetic shear flow. Eight 
periodic images of the central A^-particle system are shown. 




Figure 2: Two-dimensional version of periodic inhomogeneous boundary-driven shear flow. 
The system consists of four separate chambers of N particles each. The chambers indicated 
by arrows are driven to the right and left by moving tether forces. Heat is extracted from the 
driven chambers to maintain constant internal energy there. The boundary-driven motion of 
the other two chambers is purely Newtonian. 

Evans objected to this conclusion^, stating that the Sllod algorithm is "exact". He is 
of course correct, in the sense that the Sllod algorithm is nothing more than Newton's 
equations of motion written in a different coordinate frame, a Lagrangian frame moving 
along with the sheared fluid. But because Newton's equations by themselves cannot lead 
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to noncquilibrium steady states, the SUod algorithm is "exact'' in the rather hmited case 
of isolated systems. 

The misunderstanding evident in Evans' remark led us to undertake the present work. 
Instead of considering steady states, which seem to us the simplest situation, Evans had in 
mind a time-dependent spatially-periodic adiabatic deformation. When no thermostat or 
ergostat forces are used in the equations of motion the shear deformation is adiabatic, with 
continuous heating. In the adiabatic case no steady state results and the Sllod equations 
of motion are equivalent to Newton's equations of motion for a system undergoing periodic 
deformation with strainrate e. Just as in the thermostated case the normal stresses and 
temperatures differ. Not only the magnitudes, but also the orderings, of these components 
can, and do, differ from those found in steady states. 

In this paper we motivate and describe large-scale adiabatic-shear simulations and 
discuss the interpretation of these simulations. These simulations use periodic boundary 
conditions, just as shown in Fig. 1, but are extended here to three Cartesian space 
dimensions {x,y,z}. We take into account the important role of fluctuations in defining 
local values of the velocity, and the temperature and stress tensors. Sec. II describes the 
algorithms, and Sec. Ill the various definitions of temperature for nonequilibrium (as well 
as equilibrium) systems. Sec. IV outlines the results of the current simulations. Sec. V 
gives the conclusions we have reached as a result of this work. Sec. VI suggest extensions 
of this work. 

II. SHEAR FLOW ALGORITHMS 

Two numerical algorithms, "Sllod" and "Doll's", for spatially-periodic shear flow in 
a volume V both satisfy the macroscopic energy-balance relation E — —ePxyV. The 
corresponding solutions differ in effects of order e^, with P^^ > Pyy in the Sllod case, and 
the reverse using DolFs algorithm. To describe noncquilibrium situations it is natural to 
introduce the gradients and time derivatives of these variables, with the simplest situations 
those "stationary states" (necessarily driven by external forces or heat sources) in which all 
the partial time derivatives (the rates of change at a flxed location) vanish. Steady shear 
flow can be simulated with homogeneous sources and sinks of momentum and energy 
through the Doll's and Sllod algorithms. The adiabatic versions of these equations of 
motion (no thermostats or ergostats) introduce an overall flowfield imposed with the 
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parameter e through periodic boundary conditions: 

X = (px/m) + ey ; y — (py/m) ; [SUod or Doll's] . 

The SUod algorithm is simply a rewriting of Newton's equations of motion for the evolution 
of the "momentum" p^ relative to the motion induced by the periodic boundary conditions: 

X = (Fx/m) ; y = (Fy/m) < — ^ 

mx - mey ^ p^ ^ F^ - epy ; my ^ Py ^ Fy [SUod] . 

In the laboratory frame (where one sees the overall strainrate e induced by the periodic 
boundary conditions) the motion follows from the usual Hamiltonian, 

If, as it is in the SUod algorithm, the momentum {px,Py) is defined instead in the 
comoving frame then there is no analogous Hamiltonian. To see this in detail suppose 
that the comoving equations of motion (describing the Newtonian dynamics) could be 
derived from a hypothetical comoving Hamiltonian, T-Ccom{{x,y,Px,Py})'- 

X = +{dHcara/dpx) = (Px/m) + ey ; y = +{d'Hcom/dpy) = {py/m) ; 

Px = -{dHcom/dx) = Fx- epy ; py = -{dHcom/dy) = Fy . 

The second partial derivatives of the hypothetical Hamiltonian with respect to y and Px 
should be equal. But we find instead 

{d/dy){m,ora/dpx) = {d/dy){ey) = e , 

and 

{d/dpx){m^Jdy) = {d/dpx){-Fy) = , 

showing that there is no such comoving Hamiltonian. 

On the other hand, the very similar Doll's- Tensor equations of motion (which are not 
Newtonian) do follow from a special Hamiltonian appropriate to the comoving frame: 

^Doll's = ^p^/{2m) + ^ + e^ypx : 

+ {d'Hi)o\vs/dpx) = x = (px/m) + iy ; +{dnDon's/dpy) = y = {Py/m) . 
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-{&Hj)o\\'s/dx) = = ; -{dHBoii's/dy) = py = Fy - ep^ . 

Both the foregoing Sllod and Doll's sets of motion equations are adiabatic, so that 
the systems they describe heat due to viscous shear as time goes on. Additional time- 
reversible frictional forces of the form —(p can be added to either set of motion equations 
to keep the energy or the temperature constant^'-: 

{ AF = -Cp ; Px = m{x - ey) ; Py = my ; p^ = mz } . 

The frictional forces make it possible to explore a spatially-homogeneous nonequilibrium 
steady state with definite values of the (time-averaged) stress and temperature. Both 
these nonequilibrium properties need proper definitions. We consider several alternative 
definitions of temperature in the following section. 

III. DEFINITIONS OF TEMPERATURE 

In statistical mechanics a longstanding definition of temperature has been kinetic, 
based on the physical picture of an ideal-gas thermometer-i^. Measuring the momenta {p} 
relative to the comoving frame of the kinetic thermometer, the kinetic-theory definition 
is: 

kTxx = {pl/m) ; kTyy = {pl/m) ; kT,, = {pl/m) . 

A simple mechanical model capable of measuring all three temperatures simultaneously 
is a dilute gas of parallel hard cubes^. 

There is also a configurational analogii^i^, 

kTxx = {F^)/{Vln) ; kTyy = {F;)/{Vln) ; fcT,, = (Ff )/(V^H) . 

The configurational temperature has no clear connection to a physical model of a ther- 
mometer, but follows instead^i from a formal integration by parts of the canonical average 

of v^n. 

(V^H) = {{Vnf/kT) . 

Both temperature definitions, kinetic and configurational, have associated ambiguities: 
[1] fluctuations in the kinetic case, and [2] rotation in the configurational case. Consider 
fluctuations first. The local velocity fluctuates in time. The thermal momentum in the 
kinetic definition has to be measured in a "comoving" frame. Once the velocity is a local 
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quantity, as well as a time-dependent quantity, its definition becomes crucial. Here we 
adopt a modification of the "smooth-particle" definition of local velocity^: 



v[r) 



w{\r — ri\ < h)vi/ w(|r — rj| < h) 



where Vi is the velocity of Particle i, and that particle lies within the range h of the 
smooth-particle weighting function w{r < h). 

SPAM [Smooth Particle Applied Mechanicsi^] provides spatially very smooth mate- 
rial properties (such as density, velocity, stress, and energy) with two continuous spatial 
derivatives. The definitions of these properties require a smooth weighting function, 
w{r < h), which must be continuously twice differentiable, normalized, and which must 
also have a finite range h. Here we adopt the simplest such weighting function meeting 
these requirements, Lucy's. In three dimensions Lucy's form for w is 

105 

WLucy = ,o [l - 6a;^ + 8x^ - Sx^] ; x = \r\/h ; 



Airr w{r < h)dr = / Anr w{r < h)dr = 1 . 
Jo 

In the following section we show that the smooth particulate velocity fluctuations mea- 
sured as temperature are best defined through a slight modification of the smooth-particle 
approach, in which the "self" contributions to the particle sums, ^ w and ^ wv, are ab- 
sent. This modification reduces the number-dependence inherent in comparing atomistic 
simulations to continuum predictions. 

At first sight, the configurational definition of temperature has an advantage over the 
kinetic one in that a calculation of the stream velocity is not required. But the current 
work led us to recognize a difficulty in defining configurational temperature away from 
equilibrium. Consider rotation. Particularly in turbulent flows, rotation is important. 
Although configurational temperature has been touted as a way to avoid defining a local 
velocity—, it also contains a small and subtle ambiguity — configurational temperature 
depends on rotation rate. 

A rotating rigid body generates centrifugal forces of order uj'^r (offset by tensile forces) 
at a distance r from the center of mass. The tensile forces contribute to the configurational 
temperature definition, 

kTc = {F^)/{V^n) , 

while the centrifugal ones do not, so that perimeter particles are apparently "hotter" than 
the cooler interior by (relatively-small) contributions of order u^. 



To sec this in a simple two-dimensional example, consider the point (x, y) viewed from 
an (X, Y) coordinate system rotating counterclockwise at the angular frequency uj: 

X — xcos{u}t) + ysm.{ujt) ; Y = ycos{ut) — xsm{ut) . 

Two time differentiations, evaluated at time t — 0, give the Coriohs, centrifugal (rotating 
frame), and centripetal (laboratory frame) forces: 

X + 2ujy - uP'x = X + "luY + up-X ; 

Y = y- 2ux -u^y = y- 2uX + uj'^Y . 

For rigid rotation at an angular velocity uj a particle at (x, y) — (r, 0) with laboratory 
frame velocity (0, a;r) and acceleration {Fx,0)/m = (— a;^r, 0) the X equation becomes: 

X ^x + 2ujy - up-x = (F^/m) + 2u?r - u;^r = = {F^/m) + uj'^r , 

showing that the atomistic force {Fx/m) exactly offsets the centrifugal force uj'^r. Thus 
the configurational temperature for rigid rotation is proportional to r'^uj^. 

The numerical work described in the following section, for relatively gentle shear flows, 
supports the view that kinetic temperature is both simpler and better behaved than 
conflgurational temperature, with smaller fluctuations in both space and time. In a work 
still in progress we contrast the two approaches for the problem of a strong shockwave. 

IV. SIMULATIONS AND RESULTS 

Throughout, we focus on a dense fluid in three space dimensions, with the short-ranged 
"soft-sphere" potential of Refs. 1 and 2 : 

0(r < 1) = 100(1 - T^Y ; * = E '^(l^^^l) ' 

where the sum over pairs includes all particle pairs within a distance unity. The total 
energy of the system consists of a kinetic part in addition to the potential energy 

E = K + ^- K = ^p7(2m) . 

We focus on the dense-fluid state of Ref. 1, with a density and energy per particle of 
unity, 

E/N = Nm/V = N/V = 1 ; 6^ = 216 < < 2744 = 14^ . 
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0<time<12 0<time<12 

Figure 3: Overall adiabatic temperature variations for adiabatic shear flows with 64 x 64 x 64 
soft spheres with an initial kinetic temperature of 0.01. The strainrate dux/dy = e is 0.5. 
{Txx,Tyy,Tzz} are plotted here. For times greater than 3 neither algorithm shows significant 
differences between the temperatures on the scale of the plots. 




Sllod 



0<time<12 0<time<12 

Figure 4: Overall adiabatic pressure variation for adiabatic shear flows with 64 x 64 x 64 
soft spheres with an initial kinetic temperature of 0.01. The strainrate dux/dy = e is 0.5. 
{Pxx^Pyy^Pzz} are plotted here. For times greater than 3 neither algorithm shows significant 
differences between the pressures on the scale of the plots. 



Data for homogeneous isoenergetic Doll's and Sllod simulations are given in Ref. 1 
along with complementary results for boundary-thermostated flows. That study showed 
that the normal stress differences in the homogeneous simulations are very different to 
those found in boundary- driven flows. The number-dependence in the temperatures and 
normal stress differences of the homogeneous flows is nearly negligible, no more than 
once the number of particles is a few hundred. By contrast, the boundary-driven 
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temperatures are quite different, as the midstream temperature increases as N'^^^. 

Here we consider in addition adiabatic deformation, with the Newtonian motion driven 
by shearing boundary conditions and without any thermostat or ergostat forces. The 
initial state is a cubic lattice with a kinetic temperature of kT = 0.01 and an initial 
energy per particle of E = K + ^ = O.OISA^ + (because the nearest-neighbor separation 
is initially unity, just beyond the range of the repulsive forces). We compute and compare 
two different kinetic temperatures, each with the three components {Txx,Tyy,Tzz}- The 
time- averaged temperature, kT^^ is 

kT^J^ ^ {mf) ; 
kTj,"^ ^ {mz') , 
while the inst ant ant aneous temperature, kT^'^^^ is 

kT-:'^{m{x,-Vx{r„t)f) ; 

kT^^'^{m{m-Vy{r„t)f) ■ 

kTlT^{m{zi-Vz{r„t)y) , 

where the instantaneous velocity at Particle i's location is a modified version of the usual 
smooth-particle average^: 

v{ri,t) = ^w{r,j)vj/ ^w{rij) . 

For simplicity we choose Lucy's weight function with the range h = 3 : 

105 

w(r < h) = -^(1 — 6x'^ + Sx"^ — Sa;"^) ; x = r/h , 

for the evaluation of all the smooth-particle sums. 

In smooth-particle simulations^i the "self terms", in the two particle sums, WuVi and 
Wii are always included. In analyzing equilibrium molecular dynamics simulations for 
local velocity fluctuations (the usual kinetic temperature) numerical work shows that 
there is a much better correspondence with equilibrium temperature when the self terms 
are omitted. We have followed that practice here. If the self terms are included in 
computing the local stream velocity the resulting kinetic temperatures are roughly 10% 
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9 < time < 10 9 < time < 10 

Figure 5: {Txx,Tyy,Tzz} are plotted here for portions of 64,000-particle adiabatic simulations 
of Figs. 3 and 4. The heaviest lines show the laboratory-frame kinetic temperature; the medium 
lines show kinetic temperature relative to the instantaneous smooth-particle velocity. The light 
dashed lines show the configurational temperatures, which fluctuate more wildly than the ki- 
netic temperatures. In the steady-state simulations of Ref. 1 the Doll's kinetic temperatures 
{0.496, 0.508, 0.493} and the Shod kinetic temperatures {0.507, 0.497, 0.493} correspond to an 
internal per-particle energy of exactly unity. 




9 < time < 10 9 < time < 10 

Figure 6: {Pxx, Pyy, Pzz} are plotted here for portions of the adiabatic simulations of Figs. 
3 and 4 (using laboratory-frame kinetic contributions) and correspond to the heavy, medium, 
and light lines respectively. In the steady-state simulations of Ref. 1 the Doll's pressures (using 
laboratory-frame kinetic parts) {2.496, 2.528, 2.482} and the Shod pressures {2.516, 2.509, 2.484} 
correspond to an internal pcr-particle energy of exactly unity and an average temperature of 
about 0.5. The shear stress is about the same for the two algorithms, Cxy — 0.343. 
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lower. With the self terms omitted the three temperature definitions, time-averaged 
kinetic, instantaneous kinetic, and configurational, all give similar results. 

The need for excluding the "self" contributions can be rationalized by considering 
an equilibrium particle i at location with velocity Vi. With its neighbors' velocities 
uncorrelated (as they are at equilibrium) the smooth-particle velocity at n is, on average. 



so that the temperature, based on the velocity fluctuations as measured by the differences, 
{vi — {vspAMifi})} , is reduced by a factor of [1 — w{0)]'^. Accordingly, we have excluded 
the "self terms" in the kinetic parts of the temperatures and pressures illustrated in the 
figures. 

Figs. 3 and 4 show the overall increase of temperature and pressure beginning with a 
homogeneous cubic crystal, at a kinetic temperature of 0.01, and ending at a homogeneous 
shearing fiuid state with a temperature somewhat greater than 0.5. The details of the 
temperature and pressure for 64,000 particles, in the vicinity of kT ~ 0.5, are shown in 
Figs. 5 and 6. The fiuctuations in the data can be reduced by using even larger systems. 
Compare Figs. 5 and 6 with corresponding results for 262,144 particles, shown in Figs. 
7 and 8. In these latter simulations the internal energy per particle reaches unity for the 
Sllod algorithm at a time of 9.799, and for the Doll's algorithm at a time of 9.480. 

At a fixed strainrate of 0.5, small-system fiuctuations can completely obscure the or- 
derings of {Ta} and {Pa}. By increasing the system size it is possible to verify that 
the transient fiuctuating temperatures and stresses in adiabatic deformation are close to 
those of the isoenergetic periodic shears, with the orderings y > x > z for Doll's and 
X > y > z for Sllod. Neither algorithm reproduces the boundary- driven ordering (at the 
same density, strainrate, and energy) x > z > y. We discuss this finding in the following 
Section. 

V. NONEQUILIBRIUM CONSTITUTIVE RELATIONS 

Models for continuum mechanics follow from conservation of mass, momentum, and 
energy. The differential expressions of these conservation relations are the continuity 

equation, the equation of motion (which introduces the pressure tensor P as the comoving 
momentum fiux), and the energy equation (which introduces the heat fiux vector Q as 
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9 < time < 10 9 < time < 10 

Figure 7: {Txx,Tyy,Tzz} are plotted here for portions of 262,144-particle adiabatic shear simula- 
tions. The heaviest lines show the laboratory-frame kinetic temperature; the medium lines show 
kinetic temperature relative to the instantaneous smooth-particle velocity. The light dashed 
lines show the configurational temperatures, which fluctuate more than do the kinetic tempera- 
tures. The kinetic temperatures are nearly the same as those in the stationary shear simulations 
ofRef. 1. 




9 < time < 10 9 < time < 10 

Figure 8: {Pxx, Pyy, Pzz] are plotted here for portions of 262,144-particle adiabatic shear simu- 
lations, and correspond to the heavy, medium, and light lines respectively. The configurational 
parts of the pressure are slightly, but significantly, larger than those found in the stationary 
shear simulations of Ref. 1. 

the comoving energy flux): 

p = -pV -v ] pv^-V -P ] pe = -V-y : P - V • Q . 

The time-and-spacc dependent state variables of hydrodynamics are taken from equilib- 
rium thermodynamics, extended to the case in which gradients and time dependence can 
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occur. The state variables at the location r and time t are the density, velocity, and 
energy, {p{r,t),v{r,t),e{r,t)} and it is assumed that the pressure and heat flux can be 
defined in terms of the present values, the gradients, and possibly the past histories, of 
these same variables. 

In the present work we have seen that both the temperature (extended from the scalar 
thermodynamic variable to tensor values) and the stress can differ for two systems with 
identical densities, strainrates, energies, and constitutive relations (because the underlying 
particles are the same). Evidently both temperature and stress depend upon additional 
state variables. The relative independence of the normal stress differences to the system 
sizei L suggests that the discrepancy between periodic and boundary-driven systems is 
insensitive to second derivatives, {VVp, VVf , VVe}, all of which vary as From 
the constitutive standpoint it is simplest to imagine a dependence of the normal stress 
differences and the temperature tensor on the rate of heating e. Such a dependence could 
be used to describe the deviation of the adiabatic transient flows from corresponding 
stationary flows. Finding an additional independent constitutive variable to distinguish 
stationary boundary-driven flows from stationary homogeneous flows is a challenging re- 
search goal. 

VI. SUMMARY 

Evans' emphasis on the exactness of the Sllod algorithm (restricted to adiabatic flows) 
is conflrmed here, as Sllod is nothing but Newton in a different coordinate frame. But 
it must be noted that the large-system adiabatic pressure tensor exhibits clear differ- 
ences from the stationary pressure tensor at the same energy, density, and strain rate. 
Sufficiently large systems, with hundreds of thousands of particles, show that the nonequi- 
librium temperature tensors of adiabatic transient flows are very similar to those of ho- 
mogeneous periodic stationary flows. The "realism" of the adiabatic flows is questionable 
because real boundaries, which normally drive, constrain, and cool flows, are absent. The 
diffusion time for an A^-particle system driven by a strainrate incorporated in its periodic 
boundary conditions varies as N"^^^, so that simulation results depend increasingly upon 
their initial conditions as system size increases. 

An interesting flnding of the present work is that the smooth-particle calculation of 
local velocity (needed for the computation of the local temperature), f spam = Yl Yl ^5 
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is best modified by omitting the self terms in both sums. Considering the fiuctuations 
in V — (vspam) in an equilibrium system, we see that omitting the self terms results in 
an exactly correct temperature, while including them leads to errors of the order of ten 
percent. 

Prom the constitutive standpoint it is simplest to "explain" the difference between the 
adiabatic and boundary-driven nonlinear properties through a dependence on e, where e is 
the internal energy per unit mass. Though the configurational temperature tensor avoids 
the problem of defining a local stream velocity it still includes the effect of rotational 
contributions, giving rise to "temperature gradients" based on centrifugal forces in the 
absence of heat flow. Fortunately these rotational temperature contributions are small, 
of order e^. 

VII. WHAT TO DO? 

One referee asked us Lenin's famous question with regard to this work's consequences. 
We hope to stimulate further investigations of microscopic systems from macroscopic 
points of view. The microscopic analogs of macroscopic temperature, stress, and fluctu- 
ations are imperfect, but vital in drawing macroscopic conclusions from particle simula- 
tions. There is much to do in understanding this correspondence better for more complex 
systems with rotational and vibrational degrees of freedom. Simulations and theories of 
elongational flow have led to unresolved controversies as to the "right way" to simulate 
such flows. See, for instance Refs. 28-36 cited in our Ref. 1. We believe that the non- 
linear aspects of steady deformational flows deserve more study. For unsteady flows even 
an "exact" algorithm such as SUod, depends in an essential way on the initial conditions 
unless the deformation rate is very small. 

Shockwaves provide more extreme tests of the correspondence between microscopic 
and macroscopic models. The significance of temperature for quantum systems away 
from equilibrium needs elucidation too. We are confident that progress along all of these 
lines can best be achieved by carrying out, analyzing, and comparing series of simulations 
such as those described in the present work. 
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